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We present a numerical investigation of the gravitational collapse of horizon-size density fluctua- 
tions to primordial black holes (PBHs) during the radiation-dominated phase of the Early Universe. 
The collapse dynamics of three different families of initial perturbation shapes, imposed at the time 
of horizon crossing, is computed. The perturbation threshold for black hole formation, needed for 
estimations of the cosmological PBH mass function, is found to be S c ps 0.7 rather than the generally 
employed <5 C w 1/3 if S is defined as AM/Mh, the relative excess mass within the initial horizon 
volume. In order to study the accretion onto the newly formed black holes, we use a numerical 
ON ■ scheme that allows us to follow the evolution for long times after formation of the event horizon. In 

ON ' general, small black holes (compared to the horizon mass at the onset of the collapse) give rise to a 

fluid bounce that effectively shuts off accretion onto the black hole, while large ones do not. In both 
cases, the growth of the black hole mass owing to accretion is insignificant. Furthermore, the scal- 
ing of black hole mass with distance from the formation threshold, known to occur in near-critical 
gravitational collapse, is demonstrated to apply to primordial black hole formation. 
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CN ! I. INTRODUCTION 

ON . 
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Primordial overdensities seeded, for instance, by inflation or topological defects may collapse to primordial black 
holes (PBHs) during early radiation-dominated eras if they exceed a critical threshold This particular PBH 

ON ■ formation process, which is examined in this paper, occurs when an initially super-horizon size region of order unity 
overdensity crosses into the horizon and recollapses. Among the potentially observable consequences of PBHs, should 
they be produced in cosmologically relevant numbers, are thermal effects due to the Hawking evaporation of small 
O |' PBHs (manifested in the gamma ray background or as a class of very short gamma ray bursts ||) or purely gravitational 
effects such as gravitational radiation of coalescing binary PBH systems Q or contribution of PBHs to the cosmic 
;| ■ density parameter. Upper bounds on these signatures strongly constrain the spectral index of the fluctuation power 
^ | spectrum on small scales || . 

Recently, the possibility that stellar mass PBHs constitute halo dark matter has received attention in the context 
of the MACHO/EROS microlcnsing detections ||. It has been suggested that during the cosmological QCD phase 
• i-H , transition, occurring at an epoch where the mass enclosed within the particle horizon, i?h ~ i, approximately equals 
one solar mass, PBH formation may be facilitated due to equation of state effects manifest in a reduction of the PBH 
formation threshold Q). 

Every quantitative analysis of the PBH number and mass spectrum requires knowledge of the threshold parameter, 
S c , (for the specific definition used here see below) separating perturbations that form black holes from those that do 
not, and the resulting black hole mass, Mbh, as a function of distance from the threshold. In a simplified picture of 
the formation process, where hydrodynamical effects are only accounted for in a very approximate way, the universe 
is split into a collapsing region described by a closed Friedmann-Robertson- Walker (FRW) space-time and an outer, 
flat FRW universe. For a radiation dominated universe, it can be shown that this ansatz yields <5 C « 1/3, where <5 C 
is evaluated at the time of horizon crossing ||. On dimensional grounds, the natural scale for Mbh is the horizon 
mass, Mh ~ R^, of the unperturbed FRW solution at the epoch when fluctuations enter the horizon. However, these 
estimates for S c and Mbh are valid only within the limitations of the employed model which cannot account for the 
detailed nonlinear evolution of the collapsing density perturbations. 

In order to determine S c and Mbh for various initial conditions, we performed one-dimensional, general relativistic 
simulations of the hydrodynamics of PBH formation. We studied three families of perturbation shapes chosen to 
represent generic classes of initial data, reflecting the lack of specific information about the distribution and classifi- 
cation of primordial perturbation shapes. Our numerical technique, adopted from a scheme developed by Baumgarte 
et al. |l3| , is sketched in Sectio n |ll} followed by a description of the general hydrodynamical evolution of the collapse 



and the results for S c (Section |l|) and a discussion of accretion after the PBH formation (Section IV). Defined as 



the excess mass within the horizon sphere at the onset of the collapse, we find S c ss 0.7 for all three perturbation 
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shapes. A numerical confirmation of the previously suggested power-law scaling of M^h with <5 — <5 C 0, related to 
the well-known behavior of collapsing space-times at the critical point of black hole formation [ ^o| , is presented in 
Section |y|. In this framework, the PBH mass spectrum is determined by the dimensionless coefficient K and the 
scaling exponent 7, such that 

M bh - KM h (S - 5 c y . (1) 

We provide numerical results for K and 7 for the three perturbation families. These values may be used, in principle, 
to determine PBH mass functions as outlined in || . 

To introduce our numerical approach and isolate the dependence of 5 C and Afbh on the initial perturbation shape 
from the impact of the equation of state, we restrict the discussion here to the purely radiation-dominated phase 
of the Early Universe. In a separate publication, we will investigate the change of S c before, during, and after the 
cosmological QCD phase transition. 

Two other groups have, to our knowledge, published results of numerical simulations of PBH formation in the 
radiation-dominated universe |l0|] . Our work differs from theirs with regard to the numerical technique, the choice of 
initial conditions, and the analysis of the numerical data. Wherever possible and relevant, we compare our methods 
and results with those previously published. 



II. NUMERICAL TECHNIQUE 



The dynamics of collapsing density perturbations in the Early Universe are fully described by the general relativistic 
hydrodynamical equations for a perfect fluid, the field equations, the first law of thermodynamics, and a suitable 
equation of state. We use a simple radiation dominated equation of state, P = e/3, where P is pressure and e is 
energy density, as appropriate during most eras in the Early Universe. The assumption of spherical symmetry is well 
justified for large fluctuations in a Gaussian distribution pi, reducing the problem to one spatial dimension. 

For our simulations, we have chosen the formulation of the hydrodynamical equations by Hernandez and Misner 
fL2[ as implemented by Baumgarte et al. [[l3| (we omit restating the full system of equations but instead refer to the 
equations published by Baumgarte et al. by a capital "B" followed by the respective equation number fl4(|). Based 
on the original equations by Misner and Sharp [ |l5[ , Hernandez and Misner proposed to exchange the Misner-Sharp 
time variable, t, with the outgoing null coordinate, u. The line element then reads (Eq. B27) 

ds 2 = -e 2 *du 2 - 2e*e A/2 du dA + R 2 dQ 2 , (2) 

where e* is the lapse function, A is the comoving radial coordinate, R is circumferential radius, and dVt is the solid 
angle element (cf. Eq. B2). After the transformation, the hydrodynamical equations retain the Lagrangian character 
of the Misner-Sharp equations but avoid crossing into the event horizon of a black hole once it has formed. Covering 
the entire space-time outside while asymptotically approaching the event horizon, the Hernandez-Misner equations 
are perfectly suited to follow the evolution of a black hole for long times after its initial formation without encountering 
coordinate singularities. This allowed us, in principle, to study the accretion onto newly formed PBHs for arbitrarily 
long times (in contrast to earlier calculations [jl(J) and therefore predict final PBH masses. 

The Lagrangian form of the Hernandez-Misner equations allows the convenient tracking of the expanding outer 
regions in a comoving numerical reference frame. It also provides a simple prescription for the outer boundary 
condition, as explained below. The extremely low ratio of baryon number to energy density in the Early Universe 
requires a re-interpretation of the comoving radial coordinate, A, in Eq. (Bl) and the comoving rest mass density, po, 
in (B3). We re-define po as the number density of a conserved tracer particle with the purpose to define the comoving 
coordinate A as the tracer particle number enclosed within R. The variable e is then defined as the energy per tracer 
particle number density, such that the energy density is e — epo. In the ultra-relativistic limit e> 1, allowing us to 
replace "1 -I- e" with "e" in Eq.s B3, B4, B6, B14, and B38. This way, the Lagrangian coordinate A can be scaled to 
order unity together with all other variables, which is desirable for reasons of numerical stability. 

Given the definition of the radial grid coordinate, a suitable discretization of A must be found. Numerical accuracy 
dictates to deviate as little as possible from an equidistant grid partition lest numerical instabilities occurring on 
superhorizon scales severely constrain the grid size (see below). On the other hand, since AR ~ R~ 2 AA in a constant 
density medium, spatial resolution is concentrated near the outer grid boundary and is worst near the origin (where it 
is needed most) in case of equidistant AA. As a compromise between accuracy and resolution, we use an exponentially 
growing cell size of the form 

AA, ; = (1 + %)AA^ , (3) 
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where C is a constant and N is the total number of grid points. Based on the standard convergence tests for numerical 
resolution, we used N = 500 and C = 12 for the results reported below. 

The canonical boundary conditions (B18 and B40) are imposed at the origin, while the outer boundary is defined 
to match the exact solution of the Friedmann equations for a radiation-dominated flat universe. Hence, the pressure 
follows the analytic solution 

where t(Aj$) is the proper time of the outermost fluid clement (identified here with the FRW time coordinate, £frw) 
and Pq and r are the initial values for pressure and proper time. The time metric (lapse) functions in (Bl) and (B27) 
are fixed at 

e* - e* = 1 (5) 

at the outer boundary in order to synchronize the coordinate times t and u with the proper time of an observer 
comoving with the outermost fluid element, t(j4n), and thereby with £frw (note that B40 synchronizes u with a 
stationary observer at infinity which is a meaningless concept in an expanding space-time). 

Owing to the presence of the curvature perturbation, the space-time converges to the flat FRW solution only 
asymptotically. Therefore, the accuracy of imposing Eq. (Q) at the outer boundary is presumed to grow with the 
size of the computational domain, removing the grid boundary farther from the density perturbation. In particular, 
it is desirable to keep the boundary causally unconnected from the perturbed region for as long as achievable. The 



hydrodynamical evolution ensuing the collapse is highly dynamical for t < 100 to (Section IB), where to is the 
FRW time at the beginning of the simulation, corresponding to the light crossing time of a comoving distance of 
approximately 9i?h- Explicit numerical experiments showed good agreement of the numerical solution at the outer 
grid with the exact FRW solution for all relevant perturbation parameters if the grid reached out to i? max ^ 9Rh- 
Extending the grid to these large radii proved to be a non-trivial task for the Misner-Sharp part of the numerical 
scheme, needed to initialize the Bernandez-Misner computation, as will be outlined below. 

As the initial data is most naturally assigned on a spatial hypersurface at constant Misner-Sharp (or, equivalently, 
FRW) time t, one must transform the hydrodynamical and metric variables onto a null hypersurface in order to 
initialize the Bernandez-Misner equations. This can be done numerically in the way described by Baumgarte et al. 

first, the initial conditions are given on a t = const hypersurface and evolved using the Misner-Sharp equations. 
Simultaneously, the path of a light ray is followed from the origin to the grid boundary and the state variables on 
the path are stored. After the light ray has crossed the grid, the Misner-Sharp computation is terminated and the 
stored state values are used as initial data for the Bernandez-Misner equations. As a consequence, the Misner- 
Sharp calculation needs to be carried out until an initialization photon starting at the center reaches the outer grid 
boundary. For the aforementioned reasons, it is preferable to use a super-horizon size computational domain. Since 
the light travelling time over such a large grid is larger than the dynamical time for collapse to a PBB (~ to), a 
black hole already forms during the Misner-Sharp calculation before the initialization of the Hernandez-Misner grid 
is completed. In order to avoid a breakdown of the Misner-Sharp coordinate system inside the event horizon, the 
evolution of fluid elements is stopped artificially if curvature becomes large. More specifically, we found that the 
Misncr Sharp equations are numerically well-behaved if the inner grid boundary is defined as the innermost mass 
shell where the Misner-Sharp lapse function fulfills e* > 0.2. The inner boundary conditions are henceforth chosen 
as the frozen-in state variables on this shell. Despite this modification of the evolution equations the final results of 
the Bernandez-Misner calculation are unaffected if the initialization photon is far away from the collapsing region at 
the time of the boundary re-definition, i.e., if the black hole forms at late times during the Misner-Sharp calculation. 
In all cases reported here, this condition is satisfied. 

The second complication that arises in the Misner-Sharp coordinate system is related to supcrhorizon scales. In a 
nearly flat expanding space-time the square of the coordinate velocity, U 2 (where U = e^^dR/dt), and the ratio of 
gravitational mass to radius, 2m/ R, grow with increasing R. Both terms cancel identically in an exact FRW universe 
for all R. On scales much larger than the horizon, however, they become much greater than unity, and thus numerical 
noise can lead to significant errors in the radial metric function, T = (I + U 2 — 2m/R) 1 ^ 2 (cf. BI2). The positive 
feedback of errors in T, po (BI3), and m (B6) leads to a numerical instability of the Misner-Sharp equations on 
superhorizon scales. It can be controlled by solving for T, po, and m at each time step by iteration, and by imposing 
a very restrictive allowed density change of Apo/po < 5 x I0~ 4 per time step. 

None of the above mentioned problems exist in the Bernandez-Misner formulation by virtue of the time slicing 
along null surfaces: avoidance of the central singularity is guaranteed by the formation of an event horizon, and the 
superhorizon instability cannot occur because every grid point lies, by definition, on the horizon. Therefore, after 
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the assignment of initial data is completed, the integration of the fluid equations in Hernandez-Misner coordinates is 
numerically stable for arbitrarily long times and on arbitrarily large spatial domains. In order to achieve reasonable 
accuracy, however, the time step size must be restricted to values much smaller than the Courant-Friedrich-Levy 
(CFL) condition. The reason is most likely found in the numerical integration of the lapse function, e , which is only 
first-order accurate [ p"3[ . This problem is most important for collapse with initial conditions close to the threshold for 
black hole formation, which leads to the formation of very small black holes and gives rise to very strong space-time 
curvature. Decreasing the time step generally causes a decrease of resulting black hole mass, Mbh- At the numerical 
resolution chosen for this problem (see Eq. f3|) , the code is unable to follow the formation of black holes smaller than 
Mbh ~ 0.1 in units of the initial horizon mass. An adaptive mesh algorithm may be necessary to resolve this problem. 
In agreement with studies of critical gravitational collapse ( [ p^|Jl7| ), our experiments indicate that only the coefficient 
K in Eq. ([!]) is affected by the time step variation, while the scaling exponent 7 appears very robust. Nevertheless, 
these problems disappear rapidly with distance to the threshold such that convergence for larger PBH masses was 
attained. 

In addition to the test-bed calculations reported by Baumgarte et al. , we verified the accuracy of the code 
including the modified outer boundary conditions by simulating the evolution of a flat unperturbed universe. Using the 
time step restriction described above, the numerical results for all hydrodynamical variables differ from the analytic 
FRW solution by less than 10~ 3 . 



III. HYDRODYNAMICAL EVOLUTION OF COLLAPSING FLUCTUATIONS 



We have studied the spherically symmetric evolution of three families of curvature perturbations. Initial conditions 
are chosen to be perturbations in the energy density, e, in unperturbed Hubble flow specified at horizon crossing. The 
first family of perturbations is described by a Gaussian-shaped overdensity that asymptotically approaches the FRW 
solution at large radii, 



e{R) = e 



1 + A exp - 



R 2 



2(i? h /2) 2 



(6) 



Here, R is circumferential radius, Rh = 2io is the horizon length at the initial cosmological time to j and £o = 3/(327rfo), 
yielding Mh = (47r/3)eoi? 3 = to for the initial horizon mass of the unperturbed space-time. In the absence of 
perturbations, cirumferential radius corresponds to what is commonly referred to as proper distance in cosmology, 
r p = ar Cl where a is the scale factor of the universe and r c is the comoving cosmic distance. 

The other two families of initial conditions involve a spherical Mexican-Hat function and a sixth order polynomial. 
These functions are characterized by outer rarefaction regions that exactly compensate for the additional mass of the 
inner overdensities, so that the mass derived from the integrated density profile is equal to that of an unperturbed 
FRW solution: 



e(R) = e 



l + A 1 



R 2 
772 I 



3R 2 
2RI 



(7) 



and 



e(R) = 



R < V3Rh 
R > \/3R h 



(8) 



The amplitude A is a free parameter used to tune the initial conditions to sub- or supercriticality with respect to 
black hole formation. The shapes of all three perturbations are illustrated in Figure (Q). 
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FIG. 1. Shapes of the critical perturbations as imposed at the onset of the simulations: Gaussian (solid line, Eq. |^), 
Mexican Hat (dotted line, Eq. Q), and polynomial (dashed line, Eq. ^|). 

The relevant dimensionless threshold parameter, S c , for the purpose of evaluating the cosmological abundance of 
PBHs is the energy overdensity in uniform Hubble constant gauge averaged over a horizon volume (i.e., synchronous 
gauge with uniform Hubble flow condition). It is equivalent to the additional mass-energy inside i?h in units of Mb- 
We find similar values — 8 C = 0.67 (Mexican-Hat), 6 C = 0.70 (Gaussian), and S c — 0.71 (polynomial) — for all 
three families of initial data in our study, suggesting that the value S c ~ 0.7 yields a more accurate estimate for the 
cosmological PBH mass function than the commonly employed S c ~ 1/3. 

In cases where a black hole is formed, wc define its mass, Mbh, as the gravitational mass, m, enclosed by the 
innermost shell that conforms to e* > 10~ 10 , the temporal evolution of all shells with smaller e* being essentially 
frozen in (with regard to proper time of a distant observer). Owing to the steep rise of e* at the event horizon, 
the exact choice of the cutoff value does not affect Mbh within the accuracy reported in this work. Unless otherwise 
specified, we henceforth quote Mbh in units of the initial horizon mass, Mh, and proper time in multiples of the initial 
time, to- 




FIG. 2. Time evolution of a near-critical Mexi- 
can-Hat density perturbation with initial S — 0.6780. A 
black hole with mass Mbh = 0.37 forms in the interior. 



FIG. 3. Time evolution of a near-critical polynomial 
perturbation with initial 5 — 0.7175, forming a black hole 
with mass Mbh = 0.36. 



Figures (ph, (||), and (^) illustrate generic features of the evolution of slightly supercritical perturbations for the 
three density perturbation families, respectively. The curves display the energy density, e/eo, at constant proper time, 
r, for each mass shell (r is given in multiples of to as labeled). In Hernandez-Misner coordinates, the lack of a well 
defined global time variable corresponding to the cosmological FRW time at infinity requires this local time slicing. 
As described in |13|, we integrate dr = e*dw and store t(A,u) together with all other state variables. The curves are 
then created by plotting the energy density along the isosurfaces of r. The radial coordinate is the circumferential 
radius, scaled such that in the absence of a perturbation it may be associated with cosmic comoving radius. Further, 
the initial horizon size, i?h = 2io, is normalized to unity in the Figures. It is interesting to note that with this type of 
time slicing e(r = const) may cease to be a single-valued function of the circumferential radius, R, in cases of strong 
curvature (cf. Figure^). In particular, at the same proper time spheres containing larger baryon number (labeled by 
A) may have smaller circumferential radius than spheres containing smaller baryon number. 

In all cases shown in Figures (||),(^), and (Q), a black hole with Mbh ~ 0.37 forms. The hydrodynamical evolution 
of the three different perturbations exhibits strong similarities: initially, the central overdensity grows in amplitude 
while the outer underdensity, if present in the initial conditions, gradually widens and levels out. A black hole forms 
in the interior. Some time after the initial formation of an event horizon, material close to the PBH but outside the 
event horizon bounces and launches a compression wave traveling outward. This compression wave is connected to 
the black hole by a rarefaction region that evacuates the immediate vicinity of the black hole. The strength of the 
rarefaction differs significantly for the Gaussian perturbation shape and the mass compensated ones: while the latter 
display only a weak underdensity that quickly equilibrates, the former gives rise to a drop in energy density by three 
orders of magnitude. 

The bounce of material outside the newly formed black hole is a feature intrinsic only to black holes very close to 
the formation threshold. It effectively shuts off further accretion of material onto the newly formed PBH. As Figure 
@ demonstrates, no bounce occurs if the initial conditions are sufficiently far above the threshold. Here, a large black 
hole (Mbh = 2.75) forms whose event horizon reaches further out, encompassing regions where the pressure gradient 
is smaller, preventing pressure forces from overcoming gravitational attraction. Slightly below the PBH formation 
threshold, a bounce occurs whose strength is proportional to the initial perturbation amplitude (Figure^), indicating 
that the fluid bounce is strongest for perturbations very close to the threshold. It is likely that previous studies JlCj ] 
failed to observe this phenomenon because their initial conditions were insufficiently close to the black hole formation 
threshold. Their numerical simulations may also not have followed the hydrodynamical evolution for long enough 
after the initial formation of the PBH. 
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FIG. 4. Time evolution of a near-critical Gaus- 
sian-curve perturbation with initial S = 0.7015, forming 
a black hole with mass Mbh = 0.37. 




2 4 6 8 

R (2 To r (t/t )- /z 

FIG. 5. Time evolution of an undercritical Gaus- 
sian-curve perturbation with initial S = 0.7006. No black 
hole is formed. 
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IV. ACCRETION 



Accretion onto PBHs and their resulting growth in mass has been a highly debated subject since the suggestion 
that PBHs may grow in proportion with the cosmological horizon mass (l) . Both analytic ]l8j and previous numerical 
studies [[l0| came to the conclusion that the growth of PBH masses by ongoing accretion is negligible except, possibly, 
for very contrived initial data for the perturbations. 
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FIG. 6. Time evolution of an overcritical Gaus- 
sian-curve perturbation with initial S = 0.7196, forming 
a black hole with mass Mbh = 2.75. 
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FIG. 7. Growth rate of Mbh as a function of proper 
time immediately outside of the event horizon, providing 
a measure for the accretion rate onto the newly formed 
black hole. The lines correspond to the black holes de- 
scribed in Figures ^| (dotted line) , ^ (dashed line) , ^| (solid 
line), and M (dashed-dotted line). 



Our results generally confirm this statement for small PBHs, but we find noticeable differences between collapse 
simulations that exhibit the fluid bounce and those that do not (Figure ^) . The rarefaction following the outgoing 
density wave efficiently cuts off the flow of material into the black hole. Comparing Figures (0) and (f|), it is recognized 
that the secondary phase of mass growth for the Gaussian shape calculation may correspond to the rise of the second 
wave crest of the strongly damped density oscillation at the black hole event horizon. This second rise in density 
is absent in the weaker bounces of the Mexican-Hat and polynomial-shaped perturbation simulations. The large 
(Afbh = 2.75) black hole, on the other hand, continues to grow at a slowly decreasing rate for long times without 
gaining a considerable amount of mass in the process. Based on these results, we expect accretion to be insignificant 
for the determination of Afbh, at least for the types of perturbations investigated here. 



V. SCALING RELATIONS FOR PBH MASSES 

Choptuik's discovery pof of critical phenomena in gravitational collapse near the black hole formation threshold 
started an active and fascinating line of research in numerical and analytical general relativity (for recent reviews, see 
pl[l ). For a variety of matter models, it was found that the dynamics of near-critical collapse exhibits continuous or 
discrete self-similarity and power law scaling of the black hole mass with the offset from the critical point (Eq. |l|). 
In particular, Evans and Coleman fl6|| found self-similarity and mass scaling in numerical experiments of a collapsing 
radiation fluid. They numerically determined the scaling exponent 7 sa 0.36, followed by a linear perturbation analysis 
of the critical solution by Koike et al. ]l7| that yielded 7 « 0.3558. 

Until recently, it was believed that entering the scaling regime requires a degree of fine-tuning of the initial data 
that is unnatural for any astrophysical application. It was noted || that fine-tuning to criticality occurs naturally 
in the case of PBHs forming from a steeply declining distribution of primordial density fluctuations, as generically 
predicted by inflationary scenarios. In the radiation-dominated cosmological epoch, the only difference with the fluid 
collapse studied numerically by Evans and Coleman []l6f is the asymptotically expanding, finite density background 
space-time of a FRW universe. Assuming that self-similarity and mass scaling are consequences of an intermediate 
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asymptotic solution that is independent of the asymptotic boundary conditions, Eq. (|lj) is applicable to PBH masses, 
allowing the derivation of a universal PBH initial mass function Q. Furthermore, cosmological constraints based 
on evaporating PBHs are slightly modified as a consequence of the production of not only horizon-size PBHs, as 
previously assumed, but the additional production of smaller, sub-horizon mass black holes at each epoch p3]. 
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FIG. 8. Black hole masses as a function of the distance to the formation threshold, S — 5 C , for three different perturbation 
shape families. Using the smallest six black holes of each family, the best fit parameters to Eq. (|l]) are: 7 = 0.36, K = 2.85, 
5 C — 0.6745 (Mexican-Hat perturbation, triangles); 7 = 0.37, K = 2.39, S c = 0.7122 (polynomial perturbation, crosses); 
7 = 0.34, K = 11.9, S c = 0.7015 (Gaussian-curve perturbation, diamonds). 

Figure (^) presents numerical evidence that mass scaling according to Eq. ([!]) occurs in the collapse of near-critical 
black holes in an asymptotic FRW space-time, and therefore applies to PBH formation. All three perturbation 
families give rise to scaling solutions with a scaling exponent 7 ss 0.36. Only the smallest six black holes of all families 
were included to obtain the numerical best fit quoted in the figure captions. On larger mass scales, deviations from 
mass scaling with a fixed exponent become noticeable; in all cases, 7 tends to increase slightly for larger Mbh- Owing 
to resolution limitations discussed in Section ([n]), we were unable to compute the formation of smaller black holes 
than the ones shown in Figure (||). 

To linear order, the scaling relation ([!]) is invariant under transformations of the control parameter 5 up to a change 
of the coefficient K. This was tested explicitly by choosing different definitions of 5 (the perturbation amplitude A, 
the total excess mass for the Gaussian-shaped perturbation, and the excess mass within the horizon volume) for the 
numerical fit and obtaining identical values for 7. 



VI. CONCLUSIONS 



In the general framework of primordial black hole (PBH) formation from horizon-size, pre-existing density pertur- 
bations, we numerically solved the spherically symmetric general relativistic hydrodynamical equations in order to 
study the collapse of radiation fluid overdensities in an expanding Friedmann-Robertson- Walker (FRW) universe. 
The algorithm is adopted from an implementation of the Hernandez-Misner coordinates jl2| by Baumgarte et al. 
p3[ . It allows the convenient computation of black hole formation and superhorizon scale dynamics by virtue of its 
time coordinate, chosen to be constant along outgoing null surfaces. 

One of the parameters entering the statistical analysis of cosmological consequences and constraints due to the 
possible abundant production of PBHs is the threshold parameter, 5 C , corresponding to the amplitude of the smallest 
perturbations that still collapse to a black hole. It generally depends on the specific perturbation shape at the time of 
horizon crossing. We studied three generic families of energy density perturbations, one with a finite total excess mass 
with respect to the unperturbed FRW solution and two mass compensated ones. Defining the control parameter, 5, 
as the total excess gravitational mass of the perturbed space-time with respect to the unperturbed FRW background 
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enclosed in the initial horizon volume, our calculations yield a similar threshold value for all three fluctuation shape 
families, S c w 0.7. 

We investigated features of collapsing space-times very close to the threshold of black hole formation embedded in 
an expanding FRW solution. If the initial perturbation is smaller than S c , it grows until pressure forces at the origin 
cause the fluid to bounce, creating an outgoing pressure wave followed by a rarefaction, but no black hole. Initial 
conditions slightly exceeding the threshold, on the other hand, lead to the formation of a very small black hole at the 
origin; however, the pressure gradient immediately outside of the event horizon is still sufficiently steep to force the 
fluid to bounce. The launch of a compression wave can be observed in simulations of all three perturbation shapes. 
It is strongest in the case of a pure initial overdensity, parameterized here as a Gaussian curve, where the density 
behind the pressure wave drops by three orders of magnitude. Increasing 5 to values significantly above 5 C , the bounce 
becomes weaker and finally disappears, signaling the failure of the pressure gradient at the event horizon to overcome 
gravitational attraction. 

This behavior has important consequences for the accretion onto PBHs immediately after their formation. If a 
bounce occurs, the inner rarefaction shuts off accretion almost completely before any significant amount of material 
has been accreted. On the other hand, black holes that form from sufficiently large overdensities, where a bounce 
is suppressed, may accrete at a slowly decreasing rate for a long time. Since most PBHs created from collapsing 
primordial density fluctuations with a steeply declining amplitude distribution form very close to S c Q , we conclude 
that accretion is unimportant for the estimation of PBH masses. This is in agreement with previous studies [ pj , 
albeit for different reasons. 

Finally, the previously suggested || scaling relation between Afbh and 5 — 5 C , based on the analogy with critical 
phenomena observed in near-critical black hole collapse in asymptotically non-expanding space-times poj ] , was con- 
firmed numerically for an asymptotic FRW background. For the smallest black holes in our investigation, the scaling 
exponent is 7 ~ 0.36, which is identical to the non-expanding numerical and analytical result ]l6| within our numerical 
accuracy. The parameter K of Eq. ([!]), needed to evaluate the PBH initial mass function derived in was found to 
range from K w 2.4 to K w 12. 
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Olinto, and V. Katalinic for helpful discussions. Part of this research was supported by an Enrico-Fermi-Fellowship. 
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